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Abstract. We study the two-dimensional structure of thin transonic accretion discs in the vicinity of black holes. We use 
the hydrodynamical version of the Grad-Shafranov equation and focus on the region inside the marginally stable orbit (MSO), 
r < r m . We show that all components of the dynamical force in the disc become significant near the sonic surface and especially 
important in the supersonic region. Under certain conditions, the disc structure is shown to be far from radial, and we review the 
affected disc properties, in particular the role of the critical condition at the sonic surface: it determines neither the accretion rate 
nor the angular momentum in the accretion disc. Finally, we present a simple model explaining quasi-periodical oscillations 
that have been observed in the infrared and X-ray radiation of the Galactic Centre. 

Key words, hydrodynamics - accretion, accretion disks - black hole physics 



1. Introduction 

The investigation of accretion flows near black holes (BHs) is 
undoubtedly of great astrophysical interest. Substantial energy 
is released near BHs, and general relativity effects, attributable 
to strong gravitational fields, show up there. Depending on ex- 
ternal conditions, both quasi-spherical and disc accretion flows 
can be realized. 

The studies of thin disc accretion span more than 
three de cades. Many result s were i ncluded in text - 
books JShapiro & Teukoiskvl 1 19831 llipunovl 1 19921) . 



Minor friction between rotating gas layers leads to the loss 
of angular momentum, so that the accreting matter gradu- 
ally approaches the compact object; this motion only slightly 
disturbs the vertical balance 1 and is commonly considered a 
small perturbation. To account for this phenomenon, standard 
model introduces a phenomenological proportionality coeffi- 
cient ass < 1 between the viscous stress tensor £ and the pres- 
sure P, = assP. As a result, the radial velocity v r can be 
evaluated as 



Lvnden-Bell ( 1969) was the first to point out that super- 
massive BHs surrounded by accretion discs could exist in 
galactic nuclei. Subsequently, a theory for such discs was 
developed that i s now ca l led standard model, or the model 
of the g-disc iShakural 1 19721 IShakura & Sunvaevl Il973t 
iNovikov & Thornell973ir 

According to standard model, if the accreting gas tempera- 
ture is low enough for the speed of sound c s to be much lower 
than the Keplerian rotation velocity vk, matter forms a thin 
equilibrium disc and moves in nearly circular orbits with the 
Keplerian velocity Vk(/) = (GM/r) 1 ^ 2 , where M is the black 
hole mass and G is the gravitational constant. The disc thick- 
ness within this model is determined by the balance of gravita- 
tional and accreting matter pressure forces, 



VK V| 



(2) 



so that for c s <K vk the radial velocity is much lower than both 
the Keplerian velocity and the speed of sound. 

The foregoing is valid far from the compact object where 
the relativistic effects are unimportant. As for the inner regions 
of the accretion disc, the general relativity effects lead to at 
least two new qualitative phenomena. First, there are no stable 
circular orbits at small radii, r < r ms for a non-rotating BH; 
r ms = 3r g is the radius of the marginally stable orbit (MSO) 
and r g - IGM/c 1 is the BH gravitational radius. This means 
that the accreting matter, that falls into the region r < r ms , 



vk 



(1) 



1 Let the vertical, or transversal, direction be a 0-direction in the 
spherical system of coordinates with the half-line 6 = being aligned 
with the angular momentum vector of the infalling matter. Then, the 
longitudinal direction is a radial one. 
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[ llgumenshchev & Belo borodovl l 1997: Balbus & Hawlev 
i lKrolik & Hawlevl2002l) .~ 



rapidly approaches the BH horizon, 2 no matter if the viscos- 
ity is present. Second, thin disc accretion on to BHs is always 
transonic. This results from the fact that, according to (0, the 
flow is subsonic outside the MSO, r > r„,„ while at the horizon 
r — r g the flow is to be supersonic (Beskinl [l997l) . Therefore, 
the study of inner accretion disc regions requires consistent 
analysis of transonic flows. 

Transonic accretion on to BHs has been the subject of 
many research pub lications. Of crucial impor tance was the pa- 
per bv lPaczvnski & Bisnovatvi-Koganl(ll98ll) who formulated 
equations of motion in terms of radial derivatives only: they 
used quantities averaged over the disc thickness. Subsequently, 
many authors considered only such one-di mensional models. 
In mo st cases, the model potential of Paczvhski & Wiita 
(1980), ip g = -GMj(r — r g ), which simulates the properties of 
the str ong gravitational field of a centrally symmetric BH, was 
used ( Abramowicz et al. 1988; Papaloizou & Szuszkiewicz 
1994[ IChen et alJll997t iNaravan et alJll997t lArtemova et alJ 
2001). The Schwarzschild and Kerr metrics, that describe the 
gravitational field of axi s ymmetric BH s , were far more rare 
(Riff ert &Heroidl Il995t IPeitz & Anpll Il997t iBeloborodovl 
Il998l iGammie & Pophaml Il998albl As for the two- 
dimensional structure, it was inv estigated mostly only nume- 
rically and only for thick discs JPapaloizou & Szuszki ewiczl 
[1994 
Il99f 

Therefore, for thin accretion discs basic results were ob- 
tained in terms of the one-dimensional approach. However, 
the procedure for averaging over the disc thickness, which is 
likely to be valid in the region of stable orbits, in fact re- 
quires a more serious analysis. In our view, it is the assump- 
tion that the transverse velocity vg may always be neglected in 
thin accretion discs, i.e. that the disc thickness is always de- 
termined by 0, that is the most debatable. This assumption 
is widely used, explicitly or implicitly, virtually in all papers 
devo ted to thin accretio n discs (see the papers cited above and 
lAbramowicz & Zurekll98lHChakrabartill996l) . 

lAbramowicz et all \ 19971) pointed out the importance of al- 
lowing for the transverse velocity in the inner accretion disc 
regions. They argued that the transverse component of the dy- 
namical force has the same asymptotic near the equatorial plane 
as the transversal components of the gravitational and pressure 
forces (under the natural assumption vg oc cos 6 in the limit 
cos 6 — > 0). Therefore, the fact that the transversal velocity 
vanishes on the equator by no means implies that [(vV)v] fl may 
be ignored. However, in the end the authors concluded that for 
thin discs the dynamical term may still be disregarded so that 
the accretion disc thickness can be determined from the bal- 
ance of the gravitational force and the pressure gradient up to 
the BH horizon. We show that this conclusion is inapplicable 
for the class of problems we consider below. 

To refine the issue, it may be fruitful to recall the simplest 
case of tra nsonic accret ion, namely ideal spherically symmetric 
accretion (Bondi 1952). Far outside the sonic surface in such a 
flow the dynamical term is negligible in comparison with the 
pressure gradient and the gravitational force. However, all the 



terms become of the same order of magnitude near the sonic 
surface. Moreover, after passing the sonic surface, the motion 
of matter differs little from free-fall, therefore the pressure gra- 
dient no longer affects the dynamics of accreting matter. This 
makes it seem unreasonable to determine the disc thickness 
from 0, at least in the supersonic region. 

However, equation Q can be easily modified to account for 
dynamical forces. Indeed, we have at the border of the disc 



dH 

dr 



(3) 



where = tt/2 — 6. Now assuming linear dependence of v@ 
on for <sc 1 , we obtain 



r dH\ „ 

v r 0. 

H dr 



(4) 



With this in hand the dynamical term can be written as 



u m i d(rv @ ) dv & 

[(vV)v] = v r — — + v & —- 
ror rov) 



v r d 



r dH \ I r dH \ d ( r dH \ 

h^H + (h^^)^0(h^^)- (5) 

Now, assuming = H/r and using q}, the ©-component of 
the Euler equation, divided by (-pH/r 2 ^j , can be written as 



rv r — 
dr 



r dH 
HlTr^ 



r dH 
~ [ ~H~dr~ 



V — c V 

H 2 9 



(6) 



In the dynamic time T d ~ [v r (r ms )/c] 1/3 r„/c. 



where the two terms in the l.h.s. correspond to the dynamical 
force components v r d(rvg)/(rdr) and vgdvg/(rd8) respectively. 
The r.h.s. of l|6) corresponds to the discrepancy between the 
centrifugal force v 2 . tan© jr and the pressure force dP/(rd&) 
whose absolute value is evaluated as P/H ~ pc^/H. This dis- 
crepancy was postulated to vanish in Q. 

Accounting for the dynamical forces, equation is a sim- 
ple qualitative analogue of second-order Grad-Shafranov equa- 
tion dl4> and allows significant disc thickness deviations from 
the value given by conventional prescription Q. For instance, 
nondimensionalizing equation (|6jl yields dr, the characteristic 
radial scale of change of H near the sonic surface r = r„, 
6r(r t ) ~ H(r t ). Consistently, by Q it follows that vg ~ v r 
which indeed indicates a rapid change of the disc thickness. 
Naturally, in the region of stable orbits in a similar way we get 
6r(r > r ms ) ~ r and vg ~ v r H/r <K v r which justifies the usage 
of conventional prescription outside the MSO. 

The existence of the small longitudinal scale 6r(r*) <K r g is 
one of the key results of this paper. Note also that since (|6jl is a 
second-order differential equation for H, there in fact exist two 
additional degrees of freedom, frequently overlooked before, 
and only one of the two is fixed by the critical condition at the 
sonic surface. This hints that the critical condition determines 
neither the accretion rate nor the angular momentum in the disc 
but only imposes a restriction on the function H(r), i.e. on the 
form of the flow, near the sonic surface. 

In this paper we purposefully consider extreme boundary 
conditions with radial inflow velocity at the MSO being far 
smaller than the speed of sound. This allows us to perform the 
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study of the problem analytically. Moreover, the results turn 
out to be almost independent of the radial inflow velocity, and 
in the end the observational predictions of our model depend 
only on the the black hole mass and its angular momentum and 
are totally independent of the inflow velocity and the speed of 
sound. 

Although our model might not be realized in Nature as is, 
it - as any other analytical model - bolsters our fundamental 
understanding of the physical process. Namely, this model of- 
fers an insight into the fundamental properties we might lose 
if we do not allow for the transversal velocity in thin accretion 
discs. This is the first step to understanding the physics of the 
critical condition at the sonic surface and its influence on the 
global disc structure. 

In the next section we introduce the reader to the hy- 
drodynamical version of the Grad-Shafranov equation in the 
Schwarzschild metric as well as derive some useful relations. 
In sections[3]and[4]we sequentially study the subsonic and the 
sonic surface regions. Section[5]presents the study of the super- 
sonic flow for the case of a spinning BH, i.e. Kerr metric, fol- 
lowed by the discussion of its applications to the explanation of 
quasi-periodical oscillations detected in the infrared and X-ray 
observations of the Galactic Centre (section^. 

2. Basic equations 

We consider thin disc accretion on to a BH in the region where 
there are no stable circular orbits. As argued above, the contri- 
bution of viscosity should no longer be significant here. Hence 
we may assume that ideal hydrodynamics approach is suitable 
well enough for describing the flow structure in this inner area 
of the accretion disc. Below, unless specifically stated, we con- 
sider the case of a non-spinning BH, i.e. use the Schwarzschild 
metric, and use a system of units with c = G — I. We measure 
radial distances in the units of M, the BH mass. 

In Boyer- Lindquist coordinates the Schwarzschild metric 
is JLandau & Lifshitsll987al) 

ds 2 = -a 2 df 2 +g ik dx i dx k , 



(7) 



where 

a 2 = 1 - 2/r, g rr = a' 2 , 



= r 2 , t 



2 2-2 

= vr = r sin I 



(8) 



Here and below indices without caps denote vector compo- 
nents with respect to the coordinate basis in 'absolute' three- 
dimensional space, and indices with caps denote their phys- 
ical components. The V* symbol always represents covariant 
derivative in 'absolute' three-dimensional space with g^ met- 
rics ©. 

We reduce our discussion to the case of axisymmetric sta- 
tionary flows. This makes it possible to introduce a stream 
function <&(>, 9) which defines the physical poloidal 4-velocity 
component up = u~ r e? + u- g e- g as 

emu* = — — — (V<J> x e^), (9) 

where n is the particle concentration in the comoving reference 
frame. The surfaces Q>(r, 0) = const determine the streamlines 
of the matter (O = corresponds to the streamline 8-0). 



For an ideal flow there are three integrals of motion con- 
served along the streamlines, namely entropy, 



S = S(<D), 

energy flux, 

E(<$) = pay, 
and the z-component of angular momentum, 
L(O) = /itju^, 



(10) 



(11) 



(12) 



where yU = (p,„ + P)/n is relativistic enthalpy; p m is internal 
energy density, P — nT is pressure. 

As a result th e relativistic Euler equation 
JFrolov&Novikovll998i 



1 



1 . 



nu VbQiUa) + V n f - fm(u0) —V a w + nny -V a a = 0, (13) 

m a 

in which indices a and b take on r and 6 values, can be rewritten 
in the form o f the Grad-Shafranov scalar equation for <£(r, 0) 
(iBeskinl 19971) . 



rV*<D 



1 V'OV^OV/ViOi M 2 V' k FV k <& 



a 2 cT 2 (V(D) 2 



~D J 



2cr 2 ?zr 2 (V<t>) 2 Z) 



+ -^—z — r \m 2 E a 2 L — \-l67T 3 nT — = 0, (14) 

a 2 xu 2 M 2 \ d<D d<D/ d<D v ; 

where the V. derivative acts on all the variables except the 
quantity M. The equation above contains three integrals £(0), 
L(<D), and 5 (O). Here 



D = -1 + 



64tt 4 



1 cj 

«p 1 - C V 



(15) 



\up-E 2 - a 2 L 2 - m 2 a 2 n 2 \ , (16) 



and the thermodynamical function M 2 is defined as 



At 



2 _ 47r ^ 



(17) 



Equation (I14> contains the only singular surface - the sonic 
surface which is defined by the condition D — 0. At this surface 
the equation changes its type from elliptical to hyperbolical. 

To close the system, we need to supply the Grad- 
Shafranov equation with the relativistic Bernoulli equation, 



w 2 . = y 2 — ipi — 1, which now becomes 



E 2 - a 2 L 2 /m 2 - a 2 /j 2 
a 2 fi 2 



After using definition (|9} we reach 



E 2 = a 2 n 2 + —V- + 



M 4 



m 2 64n 4 m : 



(VO) 2 



(18) 



(19) 
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For the sake of simplification we adopt the polytropic equa- 
tion of state P - k(S )n T so that tem perature and sound velocity 
can be written as ( Shap iro" & Teukolskvl 19831) 



T = k(S)n [ 



c 2 = -k(S)n r - 1 . 



(20) 



Since the disc is thin, c s <K 1 (see Q). Therefore we can write 
yi = m p +m p W, where W = c 2 /(T— 1) is non-relativistic enthalpy 
and trip is particle mass. For an ideal gas with T = const the 
function k{S ) can be shown to have a quite definite form 



k(S) = )fcoexp[(r- 1)S], 



(21) 



which for the case of the polytropic equation of state can be 
obtained from J20i and the thermodynamical relationship dP = 
nd\i - nTdS . 

Note that equation dl9l enables us to express the quantity 
Ai and together with it all the other thermodynamical quantities 
in terms of the potential 0(r, 9) and the three motion invariants. 
Consequently, equation dl4> does contain the only unknown 
function <t>(r,0). 



3. Subsonic flow 

In this section we show that the role of the dynamical terms be- 
comes dominant with approach to the sonic surface r = r,(ff). 
The problem of passing through the sonic surface and the su- 
personic flow structure are dealt with later in this paper. Here 
we limit our discussion to the subsonic region only, where the 
poloidal velocity m. is far lower than the sound velocity c s (as 

we assume in Sec. [2 the condition w?/c 2 «: 1 holds at the 
MSO, and in this section we always neglect the terms of the 
order of w?/c 2 or higher). We can simplify Grad-Shafranov hy- 
drodynamic equation (I14> in this limit by neglecting terms pro- 
portional to D l : 



64tt 4 



a 2 m 2 



M 2 



-M 2 -V k 



dE 



, dL 

■u^E ah — 



dO 



dO 



le^nT— = 0. (22) 
d<£ 



The resulting equation describes the subsonic flow and is el- 
liptical. Hence, it contains no critical surfaces, and our prob- 
lem requires five boundary conditions: three conditions deter- 
mine the integrals of motion and the two remaining ones are 
the boundary conditions for the second-order Grad-Shafranov 
equation. 

Following Sec. ^ we assume that the a-disc theory holds 
outside the MSO. We adopt the flow velocity components, 
which this theory yields on the MSO r - r ms , 3 as the first three 
boundary conditions for our problem. For the sake of simplic- 
ity we consider the radial velocity, which is responsible for the 
inflow, to be constant at the surface r = r ms and equal to uq and 
the toroidal velocity to be exactly equal to that of a free particle 



A nearly parallel inflow with small radial velocity v r 



revolving at r — r a 



u ?( r m$< ®) = -"0, 
M ©(Vs> ®) = ® M o, 
M fi(r ms ,0) = 1/V3, 



(23) 
(24) 
(25) 



where (r ms , 0) <K \u- r (r ms , 0)| accounts for the fact that the 
flow is parallel 5 at the MSO. The angular variable = 7r/2 - 
is counted off from the equator in the vertical direction. 

For the sake of convenience, we also introduce another an- 
gular variable © ms = ms (<l>(r, 0)), a Lagrange coordinate of 
streamlines. This is a function of streamlines. For each stream- 
line, it gives the streamline's ©-coordinate at the MSO r = r ms . 
In other words, both points (r, 0) and (r ms , ms ) belong to the 
same streamline. Further, noting that ms (r ms , ®) = & and us- 
ing and Q3> . we arrive at 



dO = 2na r m;i n(r ms , ® ms )u cos ms d0 n 



(26) 



Next, we assume that the sound velocity is constant at the 
surface r — r ms ; this yields the fourth boundary condition 



c S ( r ms.®) = c o- 



(27) 



In the case of polytropic equation of state this means that both 
the temperature and the relativistic enthalpy are also constant 
at the surface r - r ms , i.e.: 

r(r ms ,0) = 7o, (28) 
fi(r ms ,&) = Mo- (29) 

Therefore boundary conditions J23i - d25t and (I27> directly 
determine the invariants £X0 ms ) and L(0 ms ), 



i(0 ms ) = M)COS0 ms , 

: V8?9 and l = udr ms )r ] 



(30) 
(31) 

V3r„ = 2 V3 



where eo = aoTo 
(see footnote0J. 

Condition E = const J30i allows us to rewrite the Grad- 
Shafranov equation in an even simpler form, 

<9 2 <D cos© d I 1 d<&\ 
~dr 2 



= -An 2 n 2 ^r — - 4n 2 n 
(i 2 dO 



2 2 

r cos 



2 & T -M. 

u d<t> 



(32) 



At r — r ms the r.h.s. of equation 1321 can be shown to describe 
the transversal balance of the pressure force and the effective 
potential, whereas the l.h.s. corresponds to the dynamical term 
(vV)v which is Mq/cq times larger t han each of the te rms in 
the r.h.s. and thus may be dropped jBeskin etalJ l2002). It is 
natural therefore to choose entropy S (<t>) from the condition of 
the transversal balance of forces at the surface r = r m „ 



. cos 2 ©„ 



dS 
d©I 



T L dL 



r 2 ii 2 d© 



(33) 



4 For a free particle revolving at r 
we have M,a(r mc ) =1/ V3, a Q = a(r ] 



■ r ms around a non-spinning BH 

) = V2/3, 7o = 7(0 = V473 



Landau & LifshitJl987 ah 
5 Note that the parallel flow is not exactly the same as the radial 
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where the value of L(0 ms ) is determined by 13 H . This yields 
the last, fifth, boundary condition 



100 



5(0 ms ) = 5(O)-- T ln(cos0 ms ), 

3c 



(34) 



whence owing to ( 12 li we have at r — r ms the standard shape of 
concentration, 



"(Vs> ©) ~ n o ex P 



6c 2 



(35) 



(to be exact, n(r ms ,&) = «o(cos 0) r ^ 3e o). Now with the help 
of i26\ . invariants J30l >. 13 II . and J34-I can be expressed in terms 
of O, therefore the Grad-Shafranov equation does contain the 
only unknown function 0(r, 0). 

The disc streamlines profile obtained from numerical com- 
putation of equation ( 1321 is shown in Fig. ^ the profiles of 
radial velocity and sound speed are shown in Fig. [5] In the 
subsonic region, r» < r < r ms = 3r g , the disc thickness 
rapidly diminishes, and at the sonic surface we have H(r*) ~ 
U()/co H(r m ^) so that we cannot neglect the dynamical force 
there JBeskin et all20 02). 

This can be readily derived from qualitative consideration 
as well. For c s <K 1, i.e. for non-relativistic temperature, we 
obtain from ( 1181 and ( 1201 . 



The quantity 



w 2 (r) 



4 - afiyr 1 - a 1 \ (6-r) 3 



(36) 



(37) 



a 2 a 2 9r 3 

is the poloidal four-velocity of a free particle having zero 
poloidal velocity at the MSO. Assuming m. = c s = c* and 
neglecting w 2 in d37l . we find the velocity of sound c» in the 
sonic point r — r», = 0: 



c» 



r+ l 



(38) 



Since entropy S remains constant along the streamlines, the 
gas concentration n„ at the sonic surface slightly differs 6 from 
the gas concentration at the MSO. In other words, the subsonic 
flow can be considered incompressible to a zeroth approxima- 
tion. 7 It is important that this conclusion holds not only in the 
equatorial plane because the additional term 1/3(0^, - 2 ) 
in d36l is also of the order of c ( 2 for the range of angles cor- 
responding to the representative disc thickness, < 0d; sc ~ cq 
(see 1351 1. Since the density remains almost constant and the 
poloidal velocity increases from uq to c* ~ cq, i.e. changes over 
several orders of magnitude for m? <k c2, the disc thickness 
H should change in the same proportion due to the continuity 
equation: 

H(r t ) * -H( rms ). (39) 
c 

6 Note that the concentration, certainly, changes from one stream- 
line to another. 

7 The same is true for spherically symmetric Bondi accretion 
( Shapir o & Teukolskvll983l) . 



Fig. 1. The disc streamlines profile (real scale) obtained from 
the numerical computation of equation ( 1321 for co = 10~ 2 and 
Mo = 10~ 5 . Solid lines correspond to the range of variables 
m?/c 2 < 0.2, where the solution to approximate equation d32l 
should not differ greatly from the solution to complete equa- 
tion J14I . Extrapolation of the solution to the sonic surface re- 
gion is shown with dashed lines. In the figure we use the radial 
variable t = (r ms - r)/(r ms uo): for the above choice of mo and 
cq, the point t = 200 corresponds to r — 2.994r g . 



Thus, purely radial motion approximation is inapplicable in the 
vicinity of the sonic surface, and both components of the dy- 
namical force become comparable to the pressure gradient near 
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Fig. 2. Profiles of the poloidal velocity m. and sound velocity 
c s at the equator for cq = 10~ 2 and uq = 10~ 5 . Solid lines 
correspond to the same range of parameters as in Fig.^ 



the sonic surface JBeskin et all 2002), 



due. 



due. 



c 2 © 



(40) 



r a® or [i 
with the position of the sonic surface r» being determined by 

(41) 

where the logarithmic factor 



r ms r * ~ ^ U Q r ms» 



2/3 



In I 21 

Uq 



2/3 



5-7, 



(42) 



and the radial logarithmic derivative of concentration near the 
sonic surface being estimated as 



rji = (r/n)(dn/dr) ~ u Q l . 



(43) 



In other words, if there appears a nonzero vertical velocity 
component, the dynamical term (vV)v cannot be neglected in 
the vertical force balance near the sonic surface. This property 
is apparently valid for arbitrary radial velocities of the flow, 
i.e. even if the transverse contraction of the disc is not so pro- 
nounced. 



4. Transonic flow 

In order to verify our conclusions we consider the flow struc- 
ture in the vicinity of the sonic surface in more detail. Since 
the smooth tr ansonic flow is an alytical at a singular point 
r _ . (-) _ () (Landau & Lifshitsl 1 987b). we can write 



(44) 



n = ra„ 1 1 + 771/1 + -r]3® + . 



ms = fl 1© + fli/J© + ^a 2 h 2 ® + ^b Q ® 3 + . . . J , (45) 

where h = (r — r*)/r*. Here we assume that all the three invari- 
ants E, L, and S are given by boundary conditions (I30i . 13 H . 
and (I34> respectively. Hence, the problem needs only one more 
boundary condition. 




/ 

sonic surface 



Fig. 3. Schematics of thin disc streamlines profile around the 
sonic point. The flow has the form of the standard nozzle. Here 
x = —h, and y = 0. 



Full stream equation (I14> may be written as 



- -V* \-L-M 2 V k <S> 



64-n 4 , d£ , dL\ , dS 

+ ^— : — — I tu E a 2 L — \-16tt\T — = 0. (46) 

a 2 m 2 M 2 \ d<5 d<D/ d<D ' 

Neglecting spatial derivatives of 11 here (this can be done in 
the vicinity of the sonic surface for cq <k 1 as is shown later) 
we get an equation with a form similar to (13 2\ . but this time 
the equation can be used in the vicinity of the sonic surface. 
Now comparing the appropriate coefficients in Bernoulli Jl 81 
and full stream equation (I46> we obtain neglecting terms of the 
order of s = Uq/Cq (see the Appendix for the exact values), 



ao 

a 2 = 

m = 



r. x (r+i)/2(r-o 
\ Co 



r + i, 

1 - a? 



Uq 



2 + 



2al 
(F+l)// 2 , 
r + 1 \ a 2 



2.25, 



„2 ' 



--(T+Dr, 2 - 



r- 1 



(47) 

(48) 
(49) 
(50) 

(51) 



where a 2 = a 2 (r») ss 2/3. As we see, coefficients i47i - i51\ 
are expressed through radial logarithmic derivative 771 (which 
we specify according to J43I as the fourth boundary condition). 
They have clear physical meaning. So, ao gives the compres- 
sion of streamlines: ao = H(r ms )/H(r t ). In agreement with ( I39> 
we have ao * co/mo- Further, a\ corresponds to the slope of 
the streamlines with respect to the equatorial plane. As a\ > 0, 
the compression of streamlines finishes somewhere before the 
sonic surface, so inside the sonic radius r < r* the streamlines 
diverge. On the other hand, because a\ <sc u Q \ for r = r* the 
divergency is still very weak. Hence in the vicinity of the sonic 
surface the flow has the form of the standard nozzle (see Fig.[5}- 
Finally, since a 2 ~ 773 ~ bo ~ u^ 2 , the transverse scale of the 
transonic region H(r„) is the same as the longitudinal one. This 
means that the transonic region is essentiall y two-dimens ional 
(a well-known fact for a nozzle, see e.g. Landau & Lifshits 
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(1987b)), and it is impossible to analyze it within the standard 
one-dimensional approximation. 

Since the transonic flow in the form of a nozzle has longi- 
tudinal and transversal scales of one order of magnitude, near 
the sonic surface we have dn * 8r ± , i.e. 



6r n * H(n). 



(52) 



Hence, for thin discs this longitudinal scale is always much 
smaller than the distance from the BH, Sr\\/r* as H(r*)/r* <K 
1 . Only by taking the transversal velocity into account do we 
retain the small longitudinal scale Sr\\ <K r g . This scale is left 
out during the standard one-dimensional approach. 

We stress that taking the dynamical force into account is 
indeed extremely important. This is because, unlike zero-order 
standard disc thickness prescription Q, the Grad-Shafranov 
equation has second order derivatives, i.e. contains two addi- 
tional degrees of freedom. This means that the critical condi- 
tion only fixes one of these degrees of freedom (e.g. imposes 
some limitations on the form of the flow) rather than deter- 
mines the angular momentum of the accreting matter. 

Connecting the sonic characteristics rji = 771 with the 
physical boundary conditions on the MSO r = r ms is rather 
difficult (for this we have to know all the expansion coeffi- 
cients in J44I and J45l ~). In particular, we cannot formulate 
the restriction on five boundary conditions fl23i - d25i . ( 1271 . 
and fl34i resulting from the critical condition on the sonic sur- 
face. Nevertheless, estimate J43i of r]\ makes sure that we know 
the parameter t]\ to a high enough accuracy. Then, according to 
( 147 \ - \5l\ . all other coefficients can be determined exactly. 

Using expansions (|44) and g5ll, one can obtain all other 
physical parameters of the transonic flow. In particular, we have 

„2 n, 



u p =c, 



2 2 

c s = c ± 



1 - 2 m h + i(T - 1) -f 2 + -(F + 1)^0 



1 



0q2 



l+(T-l)T n h+-(T- 1)^0 



1 



--(r-ixr + i)^ 



(53) 



Therefore, the sonic surface, m. = c s , has the standard parabolic 
form: 

1 



h = -(F + l)T]l® Z 



(54) 



Finally, equation (I53> enables us to estimate the spatial 
derivatives of fi. Then, it can be easily shown that the terms 
containing these derivatives in (I46i are negligible for Co « 1 
and indeed can be safely dropped. Therefore the validity of ne- 
glecting the derivatives of yU in J46l > has been justified. 

5. Supersonic flow 

Since the pressure gradient becomes insignificant in the su- 
personic region, the matter moves here along the trajectories 
of free particles. Neglecting th e VqP term in the fl-compo nent 
of relativistic E uler equation (|F rolov & Novikovl [l998). we 
have (compare to lAbramowicz et aiTl997l) 



ail? — ; 1 — h (u(p) tan 



8r 



(55) 



Here, using the conservation law of angular momentum, up can 
be easily expressed in terms of radius: — 2 y3/r. We also 
introduce dimensionless functions f(r) and g(r): 



0/(r) = ru§, 
g{r) = -aiif > 0. 



(56) 
(57) 



(58) 



Using i55\ and the definitions above, we obtain an ordinary 
differential equation for f( r) which could be solved if we knew 

df = f + 12 
dr r 2 g(r) 

From dl (St we have w? — > w 2 as r — > r g . On the other hand, 
Up at c» » co for r < r». Therefore, the following approximation 
should be valid throughout the r g < r < r* region, 

g(r) w V(aw) 2 + (ac*) 2 . 



(59) 



Equation J58l > governs the supersonic flow structure for the 
case of a non-spinning BH. To get a better match with obser- 
vations (cf. Sec. |6j, we also consider a more general case of 
a spinning BH, i.e. a Kerr BH with non-zero specific angular 
momentum a. After some calculations, equation (1581 can be 
generalized to the Kerr metric with a strikingly simple form, 



C 1L 

dr 



r 2 g (r) 



(60) 



Here eo and 1$ are the specific energy and the angular momen- 
tum of a free particle ro tating at the MSO respectively (see e.g. 
IShapiro & Teukolskvll9 83). and g (r) is a straightforward gen- 
eralization of g(r) to the Kerr case. For the Schwarzschild BH 
(a = 0, eo = V8/9, k) = 2 V3) equation J60l > reduces back 

to (ED. 

Integrating i6Q\ . we obtain 



fir) = /ctan 



(61) 



where k - ^ja 2 (\ - efy + 1^. In the equation above, 7r/2 has 
been to a good accuracy substituted for the integration constant 
arctan [/(r*)/ yfic]: for r just below r», the function / is positive 
reflecting the fact that the flow diverges; then, / = corre- 
sponds to the point where the divergency finishes, and the flow 
starts to converge. 

The results of numerical calculations are presented in 
Fig. In the supersonic region the flow performs transversal 
oscillations about the equatorial plane, their frequency inde- 
pendent of their amplitude. We see as well that the maximum 
thickness of the disc in the supersonic (and, hence, ballistic) 
region, which is controlled by the transverse component of the 
gravitational force, actually coincides with the disc thickness 
within the stable orbits region, r > r ms , where standard esti- 
mate is correct. 

Once diverged, the flow converges once again at a 'nodal' 
point closer to the BH. The radial positions of the nodes are 
given by the implicit formula f(r„) = ±oo, i.e. 



— nn, 



(62) 
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Fig. 4. The structure of a thin accretion disc (actual scale) 
for co = 1CT 2 , uo = 1CT 5 after passing the MSO r = 3r s 
(Schwarzschild case). As sufficient dissipation can take place 
in the vicinity of the first node r = r\ , we do not prolong the 
flow lines to the region r < r\ . 
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u i i i i i i 

0.001 0.01 0.1 



Co 

Fig. 5. Radial positions (in the units of M) of the nodes for 
a range of initial sound velocities. Dotted, dashed, and solid 
curves correspond to the cases a — -0.5, a — 0, and a = 0.5 
respectively. Each curve relates the radial position of a node 
to a value of the initial sound velocity. Intersection points of 
these curves with the line co = const give the the nodes' radial 
positions for that particular value of cq. 



where n = 0, 1, 2, ... is the node number; the node with n = 
corresponds to the sonic surface. In this formula the sonic ra- 
dius r, = ro can be to a good accuracy approximated by 
r ms = r ms (a) for kq « co < 1 (see d4H ). The expression for 
r m ,(a) can be found in most textbooks ( Shapiro & Teukolskv 
1983). 

Note that equations J58I and (1601 are inapplicable in the 
very vicinity of the nodes where the pressure gradient cannot 
be neglected. There may also be an additional energy release 
because the shock fronts (and, hence, extra sonic surfaces) are 



inevitably to appear there. These factors can reduce the am- 
plitude of the disc thickness oscillations. Accurate analysis of 
these factors lies beyond the scope of this paper, and we do 
not prolong the streamlines to the region r < r\ in Fig. 0] 
Nevertheless, it is quite natural to suppose that the nodes posi- 
tions are not significantly affected by the dissipative processes, 
and our simple approach may give us a good qualitative un- 
derstanding of the supersonic flow structure as long as we are 
interested in general properties of the flow. 

Figure |5] shows the positions of nodes for different values 
of co (the positions do not depend on mo for uq <k cq) and the 
BH spin parameter a. 

Travel time between the nodes, which is in some sense sim- 
ilar to the half-period of vertical epicyclic oscillations, proves 
to have weak dependence on both uq and cq. This fact provides 
a means for testing the theory via observations, and we do this 
in the following section. 

6. Applications to observations 

6.1. General ideas 

Our initial interest in researching into our model's application 
to obs ervations was kindled by the publication bv lGenzel et alJ 
(2003) of the infrared observations of the Galactic Centre (GC): 
one of the flares obviously showed a number of peaks at time 
intervals chirping with time (see their Fig. 2e), and this is ex- 
actly what we would expect within our model. In this section 
we measure physical radial distances in the units of GM/c 2 and 
time intervals in the units of GM/c 3 . 

Suppose some perturbation in the disc (a 'chunk') ap- 
proaches the MSO. We expect to observe radiation coming 
from the chunk with the period of its orbital motion, 

T ms (a)=2n(r^+a), (63) 

where a is the angular momentum per unit mass of the BH and 
r ms is the estimate of the distance from the BH to the chunk. 
After a number of rotations, the chunk reaches the MSO and 
passes through the nodal structure derived earlier (cf. Sec. [5} 
generating a flare. Each time the chunk passes through a node, 
it probably generates some additional radiation, and therefore 
the flare is likely to consist of several peaks. We believe that it 
is these peaks that were discovered in the infrared and X-ray 
obser vations of the GC dGenzel et al.ll2003t lAschenbach et alJ 
2004). 

We stress that this section is quite independent of other 
parts of this work. Qualitatively, our model's only assumptions 
are the small thickness of the disc and the divergency of the 
flow in the supersonic region (the increase of disc thickness 
with decreasing radius). The latter is guaranteed provided that 
the flow has the nozzle-like structure (see Sec.0} near the sonic 
point. 

6.2. Finding the periods 

The time interval between the detection of two subsequent 
peaks (n-th and (« - l)-th) equals the time it takes for the chunk 
to pass between the two corresponding adjacent nodes (n-th 
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and (n - l)-th), T„ r) , plus the difference in travel times to the 



.001 



0.01 



. 1 



(2) 

observer for the radiation coming from these two nodes, T„ . 



T„ = T 



(i) . 



T (2). 



(64) 



n can also be thought of as the index of observed time intervals 
(counting from one). 

The last term in J64t is important when one of the two 
nodes comes close to the BH horizon: in such a case time di- 
lation makes it significantly longer for the radiation to reach 
the observer from that node which significantly increases T®\ 
However, for outer nodes this can never happen. With this in 
hand, we can estimate the value of T\ neglecting the effect of 
In the supersonic region, the particles move along almost 
ballistic trajectories, and the orbit of a particle falling off the 
MSO is the circle inclined with respect to the equatorial plane 
(so that the centre of the orbit is the singularity of the BH) and 
perturbed by the slow radial inflow motion. This is especially 
true of outermost nodes where the matter has moderate radial 
velocities and time dilation is negligible. 

Suppose that at some moment in time the particle is located 
at the sonic point (on the equatorial plane) which is very close 
from the MSO r = r ms (see i4l\ ). Having travelled for half the 
orbital period, the particle would have reached the equatorial 
plane again, i.e. the first node, therefore 



Tt £ 1/2 T ms (a). 



(65) 



This shows that the first time interval is a little shorter than half 
the orbital period at the MSO. Therefore, we would expect the 
observations to exhibit not only the rotational frequency at the 
MSO but also its double counterpart, i.e. the frequency close to 
just twice of the rotational frequency. 

Accurate calculations of T^P - the observer's time it takes 
for the matter to travel from the (n — 1 )-th node to the n-fh one - 
yield 



dr. 

(66) 



where the coefficients of the inverse metric are g" = -I?/(p 2 A) 
and g ,ip = cog"; the definitions of S, p, A, and co can be found 

| n II A (l) 

Thor ne et alJ (11986). The resulting relationship T„ (a, cq) for 

a = -0.5 is shown in Fig.[6]wifh dashed lines. 

For definiteness and simplicity, we assume that the observer 
is located along the rotation axis of the BH. On its way to the 
observer, the radiation travels along the null geodesic that orig- 
inates at a node in the equatorial plane (e.g. r = r n , 6 = tt/2) 
and reaches the observer at infinity (r = oo, 6 = 0). Using 
these as boundary conditions for null geodesies in the Kerr met- 
ric llCarterl l968). we numerically find T ( „ 2> (a, cq). 

Figure [6] shows the dependence of observed time intervals 
on the value of the speed of sound in the disc. All intervals 
T n behave very similarly: they first decrease with increasing 
Co and then abruptly increase to infinity due to time dilation. 
This increase occurs if the innermost node in the pair of nodes 
comes close to the BH horizon (this is indicated with upward 
arrows). Although each individual time interval may depend on 
Co, the range [T m i n (a), r max (a)] of observed time intervals (see 




.001 



.01 



c 



. 1 



Fig. 6. The dependence of time intervals between the peaks in 
a flare on the speed of sound in the disc, co, for a moderately 
spinning BH (a = -0.5); we measure time intervals in the units 
of GM/c 3 . Solid curves correspond to T„ - T„ + 7^ 2) whereas 
dashed ones correspond to The uppermost solid curve cor- 
responds to the time interval T\ between the 0th and 1st peaks, 
the second solid curve from the top corresponds to the time in- 
terval T2 between the 1st and 2nd ones, etc. Even though indi- 
vidual time intervals between subsequent peaks in the flare may 
depend on the temperature in the disc (which is proportional to 
c^, see 1201 ). their minimum and maximum values remain the 
same for the range of sound velocities where there are several 
intervals observed. In the particular case of a = -0.5, illus- 
trated in the figure, we have T m i n m 32 and T m . dx m 60 with all 
other time intervals lying in between. For M - 3.1 x 10 6 Mq 
we obtain r m j n w 600 s and T mRX w 1 100 s. 



the caption to Fig. |6) is independent of co for those co where 
there are several peaks observed. With such weak dependence 
on the speed of sound in the disc, we have only two matching 
parameters: the specific spin a and the mass M of the BH. 



6.3. Matching the observations 

In the flare precursor section we associate the period T ms with 
the 22 00 ± 300 s one (group 5, cf. Table 2 in lAschenbach et all 
( 200%) and t he time interval T\ with th e period of 1 1 00 ± 1 00 s 
(group 4 in lAschenbac h et a In consistency with 

the infrared observations of the flare, the periods T\, T2, etc. 

l2003h . 



chirp with the peak number (Genzeletal. ZWi), i.e. resem 
ble the QPO structure and thus form a cumulative peak of a 
larger width shifted to higher frequencies on the flares' power 
d ensity spectra (700 + 10 s, group 3, cf. Fig. 3a and 4a 
in Aschenbach et al . (2004)). We can estimate the average fre- 
quency of this peak as l/T = 1/2 (1 /T^ + 1 /T max ). 

The results of the periods' matching procedure are shown in 
Fig-El Observational data clearly rules out high positive values 
of a (i.e. the disc orbiting in the same direction as the BH spin), 
and therefore the accretion disc around the Central Black Hole 
is likely to be counter-rotating. Moreover, the data presented 
in Fig. enable us to estimate the mass of the Central Black 
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Fig. 7. Relation of M vs. a. Dashed, dash-dotted, and solid lines 
come from matching T (700 ± 100 s), T ms (2200 + 300 s), and 
To (11 00 ± 100 s) respectively. The resulting error polygon is 
bolded. 



Hole as (3 .7+ l)x 10 6 Mq. Thi s val ue is in agreemen t with the 
studies bv lSchodel etail 120021) and lChez et aflfe003l) . 



7. Conclusions 

Allowing for the vertical velocity in the vertical force balance, 
i.e. performing the two-dimensional analysis, is extremely im- 
portant, especially near the sonic surface and in the super- 
sonic region. Two-dimensional analysis brings about the ex- 
istence of the new small radial scale, which is of the order of 
the thickness of the disc, near the sonic surface. It also shows 
that the critical condition fixes neither the accretion rate nor 
the angular momentum of the infalling matter. The nozzle-like 
shape of the flow around the sonic point leads to the oscil- 
lations of disc thickness in the supersonic region that enable 
us to come up with a straightforward interpretation of quasi- 
periodical oscillati ons observed in the radiation co ming from 
the GC jGenzel et alJ2003tlAschenbach et all2004 . This con- 
strains the mass of the Central Black Hole to (3.7+ l)x 10 6 Mq 
which is in agree ment with other studies ( Sch odel et alJ f2002: 
IChez et al]l2003l) . Our interpretation also suggests that the ac- 
cretion disc around the Galactic Centre black hole is counter- 
rotating. 
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M. J. 1997, ApJ, Appendix: Exact values of expansion coefficients 

In this section we provide exact values of expansion coeffi- 
cients d47l - ( T5TV Here we account for the derivatives of the 
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metrics and keep terms of higher orders of smallness: 
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where a, = a(r») and a a = a{r ). 



